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Abstract . 

Completely implicit, noniterative, finite-difference schemes have 
recently been developed by several authors for nonlinear, multidimensional 
systems of hyperbolic and mixed hyperbolic-parabolic partial differential 
equations. The method of Douglas and Gunn or the method of approximate 
factorization can be used to reduce the computational problem to a 
sequence of one-dimensional or alternating direction implicit (ADI) steps. 
Since the eigenvalues of partial differential equations (for example, 
the equations of compressible fluid dynamics) are often widely distrib- 
uted with large imaginary parts, A-stable integration formulas provide 
ideal time-differencing approximations. In this paper it is shown that 
if an A-stable linear multistep method is used to integrate a model two- 
dimensional hyperbolic-parabolic partial differential equation, then one 
can always construct an ADI scheme by the method of approximate factori- 
zation which is also A-stable, i.e., unconditionally stable. A more 
restrictive result is given for three spatial dimensions. Since necessary 


The main results of this paper were presented at the SIAM National 
Meeting, Madison, Wis., May 24 to 26, 1978, and section 9 was part of a 
presentation at the 751st Meeting of the American Mathematical Society, 
San Luis Obispo, California, Nov. 11 to 12, 1977. 
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and sufficient conditions for A-stability can easily be determined by 
using the theory of positive real functions, the stability analysis of the 
factored partial difference equations is reduced to a simple algebraic 


gu . ■ 
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1. Introduction . 

Alternating direction lnpliclt (ADI) methods for parabolic equations 
were originated by Douglas [10] and Peaceman and Rachford [24]. A general 
procedure for constructing ADI schemes for multidimensional parabolic 
equations and the second-order wave equation was devised by Douglas and 
Gunn [12]. An ADI method for first-order linear hyperbolic systems in 
two space dimensions was constructed by Gourley and Mitchell [15]. 

Recently, completely implicit, noniterative, finite-difference 
schemes have been developed by several authors for nonlinear, multi- 
dimensional systems of hyperbolic [1,25] and mixed hyperbolic-parabolic 
12,5,6,19,25] partial differential equations. Lindemuth and Killeen 
formulated their ADI scheme by following the Douglas, Peaceman-Rachford 
procedure, Briley and MacDonald devised their ADI algorithm by a formal 
application of the Douglas-Gunn procedure, while Beam and Warming 
constructed an ADI method by using approximate factorization. 

The linear stability analysis for the algorithms applied to systems 
of hyperbolic-parabolic equations is in a very rudimentary state. The 
primary reason is that the operators involved do not commute. In addition, 
the stability analysis of schemes for mixed hyperbolic-parabolic equations 
is generally more difficult than the analysis for either type treated 
separately. This is particularly true for schemes using more than two 
time levels. In fact, we ore aware of only one stability proof for a 
multistep ADI scheme applied to a model equation with both convection 
(hyperbolic) and diffusion (parabolic) terms [7], 

The eigenvalues associated with a mixed hyperbolic-parabolic system 
(for example, the equations of compressible fluid dynamics) are often 
widely distributed with large Imaginary parts. Since the eigenvalue 
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spectrum cannot be bounded away from the imaginary axis, A-stable linear 
uultistep integration formulas provide ideal time-differencing approxi- 
mations (For the definition of A-stable methods, see, e.g., [8,14] or 
section 2.) Although the temporal accuracy of an A-stable linear multi- 
step method (LMM) cannot exceed two [8], this is compatible with the 
accuracy achievable with typical ADI schemes. 

The purpose of this paper is to show that if one uses an A-stable 
LMM to integrate an evolutionary partial differential equation (PDE) of 
the form 


( 1 . 1 ) 


3u 

3t 


(L + L )u , 
x y 


where and are linear scalar differential operators, then one 

can always construct an ADI scheme by the method of approximate factori- 
zation which is also A-stable, i.e., unconditionally stable. Since 
necessary and sufficient conditions for the A-stability of an LMM can easily 
be determined by applying the theory of positive real functions [9], the 
stability analysis of the factored partial difference equations is 
reduced to a simple algebraic test. Our most general result is for two 
spatial dimensions with a more restrictive result for three spatial 
dimensions. We should add that the stability analysis is for simple 
linear test equations and we have not dealt with noncommuting operators. 

In section 2 we briefly review the theory of linear multistep 
methods. In section 3 we describe a method of constructing an ADI 
method by starting with a linear multistep method and then using the 
method of approximate factorization. A linear (model) test equation for 
partial differential equations is defined (section 4) and then used to 
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analyse tha stability of approximate factorisation schamaa (section 5). 

The natural extension of approximate factorisation methods to three 
spatial dimensions is discussed in section 6. In section 7 we examine 
in detail the family of A-stable linear two-step methods. To illustrate 
the notions of this paper, we write out an ADI method for the three- 
dimensional heat equation (section 8). Section 9 contains an illustration 
of a reduced stability range for an approximate factorization formulation 
which does not follow the formulation of section 3. The connection 
between chis paper and the classic paper on ADI methods by Douglas and 
Gunn [12] is discussed in section 10. The final section includes a 
summary of the approximate factorization approach described in this 
paper. 


2. Preliminaries: A review of linear multistep methods (LMM) and 

A-stability . 

A linear k-step method for integrating the first-order ordinary 
differential equation 


( 2 . 1 ) 


77 - f(u) , t > 0 , u(0) - u , 

Qt O 


is defined by 


( 2 . 2 ) 


k k ivM 

E V " 1 - m E 9 j7 ' 

J-0 J -0 


where At is the step size (t ■ nAt), the coefficients a and 6^ 

are real constants with a, i 0 and not both a. and are 

k 0 0 


zero . 


The method la said to be explicit if ■ 0 end implicit otherwise. 
Consistency end normalisation ere expressed by the relations 


(2.3a.b,c) £ ■ 0 . £ J«j * 2 


e ■ i • 


It is convenient to associate with (2.2) the polynomials 


(2.4a) 


* a 

p ( c) - ^2 “j 1 * 3 * 

j-0 


(2.4b) 


t(c) - jr 

J-0 


Consequently the I.MM (2.2) can be rewritten as 


(2.5) 


p(E)u n ■ At a (E)^ • 


where the shift operator E is defined by 


( 2 . 6 ) 


_ n n+1 
Eu " u 


Linear stability of an LMM is analyzed by applying (2.2) to the 
linear test equation 


(2.7) 


where in general 


' - 'k + “t 
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la a complex constant. 
(C to the nth power), 

( 2 . 8 ) 


If one aaaumes a aolution of the form u n 
there follows the characteristic equation 

P(C) - X At o(0 - 0 , 


C 


n 


where the polynomials p and o are defined by (2.4). The stability 
region of an LMM consists of chose values of X At for which the 
characteristic equation (2.8) satisfies the root condition, i.e., its 
roots satisfy |c^| SI and the roots of unit modulus are simple. 

An LMM is said to be A-3table if its stability region contains all 
of the left half of the complex V At plane including the imaginary axis 
(Dahlquist [8]). Since the linear test equation (2.7) has a bounded solu- 
tion if and only if Re X < 0, the notion of A-st ability is equivalent to 


(2.9) stability of ODE stability of LMM 


where ODE denotes ordinary differential equation. Dahlquist [8] proved 
that the order of accuracy of an A-stable LMM cannot exceed two and that 
an A-stable method must be implicit. In addition, he showed that the 
trapezoidal formula 


( 2 . 10 ) 


n+1 

u 


n 

u * 


At/du n+ ^ du n \ 
~2\dl + dt / 


has the smallest truncation error of all the A-stable LMMs. 

The advantage of an A-stable LMM is that the stability of the ODE 
is a sufficient condition for the unconditional stability (i.e., stability 
for an arbitrary value of At) of the LMM. In this paper we extend the 
notion of A-stability to LMM methods applied to partial differential 
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equations (PDEs). The numerical scheme is constructed so that approxi- 
mate spatial factoring into a product of one-dimensional operators 
retains the A-stable property: 

(2.11) stability of PDE stability of factored LMM . 

The precise meaning of a factored method will be clarified in the 
following section. 


3. Construction of an ADI scheme using an LMM and approximate 
factorization. 

For the development that follows, a convenient form of the LMM (2.5) 
is 


(3.1) p(E)u n - uAtp(E) jjj 


At [a (E) - up (E) ] 


du 

dt 


0(At 3 ) 


where 


u 



Henceforth we consider only A-stable LMMs and assume that (3.1) is 

A-stable. Since the order of accuracy of an A-stable LMM cannot exceed 

two, we concentrate in this paper on the second-order schemes as 

3 

indicated by the symbol 0(At ) on the right-hand side of (3.1). The 
parameter u is defined so that the operator o(E) - up(E) on the 
right-hand side of (3.1) is at least one degree lower than the operator 
p (E) on the left-hand side. This is obvious by noting that 
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4* . . 


o(E) - e k E k + 


and hence 


o(E) - ^ (E) 



Consequently, the right-hand side of (3.1) can be computed explicitly 
from known data when advancing the numerical solution from n + k - 1 
to n + k. 

Insertion of the linear PDE (1.1) into (3.1) yields 


(3.2) [I - wAt (L + L )]p(E)u n - At[o(E) - u>p(E)](L + L )u" + 0(At 3 ) 
x y x y 


Here, for simplicity, we have assumed that the linear operators and 

Ly are independent of time although the case of time-dependent and/or 
nonlinear coefficients can be handled without difficulty [3], It is 
important to note that the unknown variable to be computed is p(E)u n 
and not u n+k . This choice ensures that the approximate factorization 
does not upset either the temporal accuracy of the scheme (see below) or 
the unconditional stability (section 5). 

If the spatial operators L and L are approximated by 

x y 

appropriate difference quotients, one obtains, in general, an enormous 
linear system to solve for p(E)u n . The computational problem can be 
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reduced Co a sequence of one-dimensional (inversion) problems by an 
approximate factorization of the left-hand side of (3.2): 

(3.3) (I - u>At L )(I - uAt L )p(E)u n - 

x y 

dtlo(E) - wp(E)](L + L )u" + 0(At 3 ) . 

x y 

On comparing the left-hand sides of (3.2) and (3.3), we see that they 
differ by che cross product term 

u> At L L p(E)u n . 
x y 

But by expanding p(E)u° in a Taylor series about u n and using 
the consistency and normalization conditions (2.3a,b), there follows 

p(E)u n - At |j- + 0(At a+1 ) , a £ 1 . 

Consequently, the cross product term 

2 2,, n 2.3,. 3u n , ... ci+3. 

u; At L L p(E)u ■ u At L L r + 0(At ) 
xy x y 5t 

- 0(At 3 ) 

is a third-order term and formal accuracy of the scheme (3.2) is not 
upset by the approximate factorization (3.3). The computational sequence 
to implement the factored scheme (3.3) as an alternating direction 
sequence is not unique. Perhaps the most obvious choice is 

(3.4a) (I - wAt L )p(E)u* - At(c(E) - u»p(E)](L + L )u° , 

x x y 
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(3.4b) 


(3.4c) 



where p(E)u is a dummy temporal difference. Here v have used 
notation from the theory of LMM for ordinary differential equations and 
consequently the algorithm (3.4) looks rather unfamiliar. In section 8 
we write out an example following notation more conventional for partial 
differential equations. 

An important aspect of implementing an ADI scheme is determining 
the proper boundary values for the intermediate dummy variables such as 

A 

p(E)u . Although such considerations are outside the scope of the 
present paper, they have been discussed elsewhere (see, e.g., [13,2]). 


4. A linear test equation for partial differential equations. 

For first-order ordinary differential equations, Ec. (2.7) is 
known as the linear test equation. For a first-order evolutionary PDE, 
we define a linear test equation for two spatial dimensions as 


(4.1) 


2 2 2 

3u , 3u , 3u 3 u . , 3 u 3 u 

5 t + c l ^ + c 2 T 7' a ^2 + b 3^7 + C '^2 • 


where c^, c^i a, b and c are real constants. In order to determine 
the conditions to be imposed on these constants for which the PDE (4.1) 
has a bounded solution, we seek a solution of the form 


(4.2) 


iUjX+K y) 

u(x,y,t) ■ v(t)e 
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f i *' L— J LA 


• . r . > > mj. uni 









t 


i 



where v(t) Is the Fourier coefficient, and are the Fourier 

variables (wave numbers). The Fourier coefficient v(t) satisfies 

(4.3a) 4? - Xv , 


where 

(4.3b) A - -Kc^ + c 2<2^ “ (a^ 2 + bic 1 ic 2 + CIC 2 2 ^ * 

For the PDE (4.1) to have bounded solutions, Re A £ 0. Consequently, 
the quadratic form 

2 . 2 
atc^ + + c< 2 


must be nonnegative for arbitrary values of ic^ and tc^* which 

2 

implies that a,c St 0, b 5! 4ac. In the absence of the convective terms, 

2 

i.e., “ c 2 “ tbe inequalities a > 0, b < 4ac are the conditions 

under which (4.1) is parabolic. The convection coefficients c^ and 
c^ are arbitrary real numbers, and in the absence of diffusion, i.e., 
a « b * c * 0, the eigenvalue A is pure Imaginary. 


5. Linear stability analysis for two-dimensional unfactored and 
factored schemes. 

In this section we examine the stability of the unfactored scheme (3.2) 
and its factored counterpart (3.3) when applied to the linear test 
equation (4.1) with the mixed derivative bu X y set to zero * The linear 
operators I. and 1. are 
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3C 

t 



T 



i 







I.lf.J 1 L i t 1 1 i_. J 


L • 

X 


J, 3 2 


(5.1) 


L • 

y 


j. . r 21 

‘ c 2 3y 3y 2 


A stability analysis including tha mixed spatial derivative is considered 
in appendix A. 

Ut assume for simplicity that u n is spatially continuous and assume 
a solution of the form 


(5.2) 


n 1 (k x+K,y) 
v e 


where v n is the Fourier coefficient and are ^ ourler 

variables. In practice, the spatial derivatives are replaced hv discrete 
difference quotients; however, as indicated at the end of this section, 
the stability proof for the spatially discrete case requires only a minor 
modification of the following stability proof. We first consider the 
stability of the unt'aotored scheme (1.2). Equation (3.2) was written 
in the particular form shown as a precursor to factoring the scheme Into 
a product of one-.llroenslooal operators. However, as It stands, (1.2) Is 
simplv 

c(F.)u° *■ At o(E) (L + l. )u n , 
x y 


where, for the present analysis, 1. and l are defined bv (5.1). 
i xy 

Assuming a solution of the form (5.2), we obtain the characteristic 



I ) 





equation 



(5.3a) 


p(C) - X At o(C) - 0 . 


where 


(5.3b) 


2 2 

A - + C 2*V " ^ aK l * CK 1 * 


and c is the amplification factor defined by 


v5.4) 


n+1 n 
v - Cv 


But (5.3a) has the sane form as the characteristic equation (2.8) 
obtained when an LMM method is applied to the linear test equation (2.7) 
for ordinary differential equations. Since Re A £ 0, the resulting 
scheme is unconditionally stable since we have assumed that the original 
LMM is A-stable. We note that for an unfactored scheme with u 11 assumed 
to be spatially continuous, an equivalent way of obtaining the character- 
istic equation (5.3a) would be to apply an LMM directly to the ODE (4.3) 
satisfied by the Fourier coefficient. 

For the factored scheme (3.3) where L and L are defined bv 

x y 

(5.1), we again assume a solution of the form (5.2) and obtain the 
following characteristic equation for the amplification factor: 


(5.5a) 


where 


(X x + * 2 )At 


1 + io At A A ^ 


o(0 - 0 , 


(5.5b) 


2 2 
A ^ “ a,c l • x 2 " “ ic 2 K 2 “ CK 2 
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■ i m ul 



with 


Re , Re £ 0 . 

On comparing (5.3) and (5.5) we see that X ■ A^ + has been replaced 
by 


(5.6) 


AF 


*1 + *2 
2 2 

1 + w At A^A 2 


The product term in the denominator is a direct result of the approximate 
factorization (AF) and hence we call this "eigenvalue" A^. A sufficient 
condition for the factored scheme to be unconditionally stable Is 
Re A af £ 0 which follows directly from the following lemma: 

A + A 

(5.7) ReL, ReA SO ^ Re £ 0 


for arbitrary real a. 

PROOF. Recall for an arbitrary complex number z that if Re z £ 0, 
then Re z * £ 0. Then, by noting that 


1 + a V 2 
“ A 2 


A i + A 2 


° 2 y 2 

*1 + A 2 


1 a 

a l + a 2 + 1/X L + i/a 2 


the lemma obviously follows. Q 
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In the above anelyela we assumed that spatial derivatives were 
continuous. It remains to consider the spatially discrete case. We 
conside. th< case where the spatial derivatives in (5.1) are replaced 
by thrne-po mt central difference quotients: 


(5.8) 


fill ■ ' ' + 0(lx2) > x ‘ 3 4x • 


.2 U... - 2u. + u. , 0 

M - -ill i — ii + o(4x 2 > 

3x Ax 


with analogous expressions for the y-derivatives. The stability 


analysis proceeds as in the spatially continuous case with the exception 
that the exponential in (5.2) is replaced by 


(5.9) 


1 n iOcjjAx-Hc kAy) 
u ■ v e 


where x » j Ax, y ■ k Ay. If we make the following correspondence 


(5.10) 


2 sin(8 1 /2) 


. < 2 


2 sin(0 2 /2) 


^ cos (0^/2) , c 2 


c 2 cos(0 2 /2) . 


where 


^ ^ AX , 02= <2 A y 


between the parair.-ce s for the discrete and continuous case, then the 


eigenvalue (5. Zo) has the same form for both the continuous and the 
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discrata caaa. Recall for th« spatially continuous casa that Ra A £ 0 


for arbitrary rsal valuas of c i» c 2» ic i»*2 and hence the discretisation 
does not change the essential property for A-atabillty that Re X S 0. 
Likewise* for the factored algorithm <3. 3) with discrete spatial deriva- 
tive approximations, one obtains by using the correspondence (5.10) the 
same characteristic equation (5.5) as for the spatially continuous case 
and again the algorithm is unconditionally stable. 

In this section we considered the stability of the factored scheme 
(3.3) when applied to the linear test equation (4.1) with the mixed 
derivative bu^ set t0 z * r °. In appendix A we prove that if the mixed 
derivative is treated explicitly, then the resulting factored scheme is 
unconditionally stable. The particular explicit method considered in 
appendix A has the defect that the resulting scheme is first-order 
accurate in time for the mixed derivative. Beam and Warming [2] showed 
that it was possible to construct an unconditionally stable two-step 
scheme where the mixed derivative is second-order time accurate by using 
linear extrapolation. However, the characteristic equation for the 
amplification factor did not have the form (5.5a); and the stability 

analysis did not include the convective terms c,u and c.u . 

1 1 x 2 v 


6. Extension to three spatial dimensions. 

In three spatial dimensions, an obvious generalization of the 
two-dimensional approximate factorization scheme (3.3) is 

(I - ujA t L^)(I - wAt L^MI - wAt L z )k>(E)u n - 

(6.1) At(o(E) - u)O(E)] (L x + L + L z )u n . 
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For three spatial dimensions a linear test equation (without mixed 
derivatives) is 


2 2 2 

’>\ <* u . „ 3u 3u , . 3u _ n 3 u , „ 3 u , 3 u , 

(6.2) at + c i 3x + c 2 3 y + c 3 3z ®1 3x 2 * *2 3y 2 a 3 ^2 


where the coefficients c £ are arbitrary real numbers and a^ Si 0. 

If the approximate factorization scheme (6.1) is applied to the test 
equation (6.2), one obtains the following characteristic equation for 
the amplification factor: 


(6.3a) 


where 


p( 0 - * AF At o(;) = 0 


(6.3b) X 


*1 + A 2 + X 3 


AF 1 + U 2 At 2 (X 1 X 2 + X 2 X 3 + XX) - W 3 At 3 X x X 2 X 3 


and 


-ic.<„ 
£ £ 


£ £ 


, Re \ l £ 0 


The lemma (5.6) of the previous section does not extend to three 

dimensions for arbitrary values of A^ with Re S 0. We consider 

2 

two special cases. If A^ is pure real, i.e., A^ = -a^K^ , then the 
denominator of (6.3b) is positive and, consequently, Re A = A <! 0 
and the scheme is unconditionally stable. If A^ is pure imaginary, 
i.e., X ^ * -ic ic , then A^ has the form 


r.-f 




h 


a i 


( 6 . 4 ) 


1 i; f ; fj i i i i i; \ i f r f m. i r. 


ie i 

A AF " o 2 + i0 2 ’ 


where 


B 1 " "* C 1 K 1 + c 2 K 2 + c 3 K 3 ) * 

8 2 - - U 3 At 3 c 1 c 2 c 3 < 1 < 2 < 3 , 

2 2 

<* 2 ■ 1 - (*) At (c 1 c 2 k i k 2 + C 2 C 3 K 2 K 3 + C 3 C 1 K 3 K 1^ * 

The region of stability for a typical A-stable LMM is illustrated in 
figure 1. For given values of c^.c^c^, one can always pick wave 
numbers k^,< 2 ,k 3 such that A At as defined by (6.4; has a positive 
real part and falls in the unstable "hole" of figure 1. Consequently, 
for pure imaginary eigenvalues (a^ ■ 0), the unconditional stability of the 
three-dimensional factored algorithm for the model equation (6.2) is not 
retained. 



7. A-stable linear two-step methods . 

Since the order of an A-stable LMM cannot exceed two, one- and two- 
step schemes are of primary interest. The addition of more steps or 
time levels complicates the numerical scheme and generally increases 
computer storage requirements with no attendant increase in accuracy. 

The most general consistent two-step method [i.e., k = 2 in (2.2)] can 
In* written as 



» 
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( 7 . 1 ) 


(1 + Ou n+2 - (1 + 20u n+1 + $u n - 


(•ft' 


n+2 . n+1 

1 + (i - e + +) 


"1 

j n 

"J" 1 


where the parameters 6,£,$ are arbitrary real numbers and e Is the 
local truncation error 


(7.2) 


e n - (« - $ + 6 - -bat 2 i-iL + 0(At 3 ) 

z dt 


determined by a Taylor series expansion about t » n At. The class of 

all two-step methods that are at least second-order accurate is obtained 

2 2 

by setting the coefficient of d u/dt in (7.2) to zero, 


(7.3) 


C - 0 + 2 


in which case the local truncation error is 


(7.4) 


e" - (28 - K - | ) 4- — T + 0(AtS 

dt 


If f. ■ $ ■ 0 in (7.1), we obtain the linear one-step method (which is 
a subclass of the two-step method), 


(7.5) 


, j n+1 , n 

n+1 n , ,, du . , , du n 

u - u - At 0 - + (1 - 0) - - e 


where we have shifted the time index down by one. This scheme is 
sometimes called the 0 -method. If 0 - 1/2, we obtain the trapezoidal 
formula (2.10), which is the only socond-order-aeeuruio one-step method. 


m 0§I I mm* 


•w-t 
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Because the trapezoidal formula has the smallest truncation error 
of all A-stable LMMs [8], one might ask why we bother to consider the 
class of A-stable linear two-step schemes. Unfortunately, the trapezoidal 
formula has the property that the characteristic root t -1 as 
A At ». Consequently, when applied to stiff ordinary differential 


» 

i 

t 

i 


j 


equations, the trapezoidal formula can produce slowly decaying numerical 
oscillations. When the trapezoidal formula is applied to hyperbolic 
partial differential equations where central spatial difference approxi- 
mations are used, the resulting algorithm is neutrally stable, i.e., 
the eigenvalues of the amplification matrix have unit modulus. If the 
solution is nonsmooth by virtue of the presence of shock waves, shear 
layers, etc., the resulting solution can exhibit highly oscillatory 
errors. In both of these applications, the oscillatory errors can be 
damped by a filtering procedure [18] or by the addition of dissipative 
terms in the case of hyperbolic equations [1,2]. An alternative to 
using the trapezoidal rule with smoothing is to use a "nonsymmetric" 
A-stable scheme such as the second-order backward differentiation 
formula [(7.1) with 0 ■ 1, £ ■ 1/2, <p * 0] . Nevanlinna and Liniger [23] 
have recently suggesced contractive methods for problems with nonsmooth 
solutions or where the lack of smoothness is introduced by rapidly 
varying integration step sizes. A particular example is a method called 
the contractive Adams method [23] [(7.1) with 0 * 3/4, £ - 0, $ = -1/4]. 

An elegant and simple test for A-stability can be formulated in 
terms of positive real functions. This terminology is borrowed from the 
literature of electrical engineering (see, e.g., [16], page 409). 
Dahlquist [9] has recently extended the theory of positive real functions 
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and considered applications to stability problems arising In numerical 
analysis. By applying this theory, it is easy to show that the linear 
two-step method (7.1) is A-stable if and only if the parameters (0,1, ♦> 
satisfy the following inequalities: 

0 i ♦ + 1/2 , 

(7.6) C * -1/2 , 

t S H » • 1/2 . 

The details of the analysis for obtaining these inequalities are 
described in |3], 

In particular we are interested in the class of all A-stable methods 
that are second-order accurate as determined by the condition (7.3). 

In this case two parameters remain and the inequalities (7. h) 

defining the A-stahle schemes become 

(7.7) i S 20 - 1 , i 2 -1/2 . 

The shaded region of figure 2 shows the range of the parameters (0.0 
tor which this class of methods is A-stahle. l.inlger (20) devised a 
sufficient and "almost necessary" condition for A-stuhilitv of LMMs. 

As an application of the criterion he developed, l.inlger determined the 
constraints for A-stabtlltv on the parameters that define the family of 
all linear two-step schemes that are at least second-order accurate. 

1'he shaded domain of figure 2 reproduces the result determined hv 
liniger. t/hir parameter i.»at ion of the two-step methods Jitters from 
l.iniger but it is easv to make the proper correspondence. ) 
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In numerical algorithms for partial differential equations it is 
conventional to use n + 1 as the most advanced time level. Hence we 
multiply the linear two-step scheme (7.1) by the shift operator E * to 
obtain 


(7.8) 


E _1 p(E)u n - At E _1 o(E) ~ , 


where 


(7.9) 


E~ p (E) - (1 + OA - 57 , 


(7.10) 


E o(E) - 1 + 0A + $7 , 


where A and V are classical forward and backward difference operators 
defined by 


(7.11) 


. n n+1 n _ n n n-1 
Au ■ u -u , Vu ■ u - u 


As a notatlonal simplification we denote E p(E) by the operator A: 


Au" - E~ 1 p(E)u n - [(1 + OA - Olu n 


(7.12) 


(1 + Ou n+1 - (1 + 20u n + Oi n_1 


a n , »f2 n 
Au + £<5 u , 


where 6 is the second-central difference operator 


.2 n n+1 , n . n-1 

6u»u -2u+u 


The class of A-stable, two-step, second-order methods (see figure 2) 


is quite large. Certain subclasses offer particular advantages in regard 
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Co computer storage, numerical dissipation, etc. As an example we 
consider one special subclass that has a particularly simple computational 
form. Let (6.1) be an ADI scheme where the time differencing is the 


two-step method (7.1). In this particular case u - $ 2 /a 2 " 9 ^ 1 + 
Multiply (6.1) by E to conform to the convention that the most 
advanced time level is n + 1. The computational algorithm is simplified 
when the difference operator in the brackets on the right-hand side of 
(6.1) is the identity operator. Since from (7.9) and (7.10) 


(7.13) 


E - (o (E) - wp(E)l - 1 + 


(♦ + rh) 


this becomes the identity operator if 


(7.14) 


1 + 4 


For the class of second-order methods, (0,4,<?) are related by (7.3) and 
hence (7.14) can be written in terms of (0,4) as 


(7.15) 


0 >* (4 + 1 ) ( 4 + 1/2) . 


This curve is shown in figure 2 as a dashed lino. With the notation 
(7.12) and the condition v7.14), which reduces (7.13) to the identity 
operator, the factored scheme (b.l) becomes 


(7.16) (I - *>At L ) ( I - u'.U !.)(!- ..'At 1. ) Au = 

X V z 


where w ■ 0/(1 + 4). 


At (L U U )u 
x v z 
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8. An ADI scheme for the heat equation in three dimensions. 


To illustrate the results of the preceding sections we write out 
an ADI scheme for the linear heat equation 


( 8 . 1 ) 




A aji ? 2 u 

* ■ 82 7 + 82 7 >7 ■ 


where a. £ 0. The scheme considered will be for the class of second- 
ly 


order A-stable schemes rep r esented by the portion of the curve 
9 » (S + 1)(£ + 1/2) of figure 2 in the shaded region. When applied to 
the heat equation (8.1), the factored algorithm (7.16) is 


( 8 . 2 ) 


( ■ ■■ 5)(‘ - u4t 82 p)(* • j - u Sj au " 

(“7 


> 1\ 

B 3- 3" 3" \ n 

Atla, ~ 2 + a 2 — , + a 3 -j-lu 

x ?v 3z / 


An obvious implementation of (8.2) is 


/ 3 2 \ * ( 3 ^ 7 n 

1 - . At a lAu - At a + a 0 , + a — ^ )u 

\ 1 3x7 V 1 - 3y" 3 3 z7 


(8.3a) 1 - .>At a 


( y\ 

j 1 - .At a., — ^ lAu 


(8.3b) [1 - .At a.. 


** * 

B A vi 


/ 9“ l n ** 

(8.3c) fl - .At — -lAu « Au , 


(8.3d) 


(1 + f.)u n+1 * Au° + (1 + 2:‘.)u n - 4u U * 
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where u » 0/(1 + O - £ + 1/2 and Au and Au are dummy variables. 

2 2 

If the space derivatives 3 /3x , etc., are approximated by central 
difference quotients, then the x-, y-, z-operators on the left side of 
(9.3a,b,c) each require the solution of a tridiagonal system. There is 
a well-known and highly efficient algorithm for the solution of tridiagonal 
systems (see, e.g. [17, page 55]). The final step (8.3d) is to compute the 
solution u n+1 from known values of Au n , u 11 and u n ^ . 


9. Reduced stability boundary for an alternative formulation. 

It has been demonstrated that an approximate factorization scheme 
for the model equation (4.1) can be constructed that is unconditionally 
stable if the time differencing is based on an A-stable LMM. In this 
section we show that an alternative formulation leads to a reduced 
stability range. 

Here we con der only the hyperbolic model equation 


(9.1) 


3u 3u 3u 

3t C 1 3x C 2 3y 


0 


and the linear two-step scheme (7.1). After multiplying (7.1) by E , 
one can rewrite the two-step scheme (7.1) as ( 7 . 8 ): 


« it 

(9.2) r (1 + OA - £V]u n = At [ 1 + OA + 4>v 

d t 

where A and V are forward and backward difference operators de f ined 
by (7.11). The operator on the left-hand side of (9.2) is A = E ^p(E) 
as defined by (7.12). As formulated in section 3, the unknown variable 
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to be computed is Au 11 . In this section we alter the procedure and 
take Au 11 as the unknown variable. 

If the time derivative on the right-hand side of (9.2) is replaced 
by spatial derivatives from (9.1) there follows 


(9.3) 



An approximate factorization of the left-hand side of (9.3) yields 

<»•*> ( 1+ fn c i i)( x + it\ c 2 £) au " ■ RHS<9 - 3) ■ 


where RHS (9.3) denotes the right-hand side of (9.3). In appendix B 
we analyze the stability of the factored scheme (9.4). For the class 
of all two-step methods (9.2) that are at least second-order accurate, 
the parameters (0,i;,4>) are related by (7.3). For this class of methods, 
the parameter space (0,0 for which the factored scheme (9.4) is 
unconditionally stable is shown by the shaded region of figure 3. The 
wedge shaped region to the "right" of the dotted lines shows the parameter 
space for which the LMM scheme (7.1) with condition (7.3) is A-stable 
(this region coincides with the shaded region of figure 2). Consequently, 
a fairly large class of factored schemes that are unconditionally 
stable according to the formulation of section 3 are not unconditionally 
stable if the unknown variable is taken to be Au° rather than 
E~ 1 p(E)u n ■ Au”. We should note that the Au° and the Au° formulations 
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•re Identical for the subclass of two-step schemes where 5 ■ 0, since 
in this special case Au n * Au 11 , as is obvious from (7.12). 

The simplest computational version of (9. A) is the variant where 
$ - 0. In this case (7.3) becomes £ » 0 - 1/2 and this subclass of 
schemes was considered in [ 2 ] and [ 25 ]. The line £ ■ 9 - 1/2 is shown 
in figure 3 and falls in the region of unconditionally stable schemes for 
0 Z 1 / 2 . 

The reduced stability range illustrated in figure 3 results from 
applying the linear two-step method (9.2) to the model hyperbolic 
equation (9.1) and taking Au u to be the unknown variable. If we had 
followed exactly the same procedure for the parabolic model equation 


3u 

3t 




* 


the resulting approximate factorization scheme in the Au 11 variable 
would not have resulted in a reduced stability range in the parameter 
space (6,5). The example of this section and the result of .section 6 
for three spatial dimensions show that maintaining unconditional stability 
of approximate factorization schemes is more difficult for hvperbolic 
equations than for parabolic equations. 


10. The relation between the Douglas-Cunn method and the method of 
approximate factorization 

A general procedure for devising ADI schemes from fully implicit 
schemes for parabolic equations and the second-order wave equation was 
developed by Douglas and (.unn [12]. In this section we briefly discuss 
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Che relation between the method of Douglas and Gunn and the method of 
approximate factorization as formulated in section 3. 

Linear one-step ...ethods [defined by (7.5)] are called two-level 
difference schemes by Douglas and Gunn. For linear one-step methods, 
the operator p(E)u n is simply the forward difference operator 

(10.1) p(E)u - Au » u - u 

A method where the increment p(E)u n - Au 11 is taken to be the unknown 
variable (rather than u° +1 ) is sometimes said to be in the "delta" form. 
For a linear one-step method, the time-differencing method for the ADI 
scheme (8.3) with £ « 0 corresponds to the trapezoidal formula (2.10). 

In this special case, the scheme (8.3) for the heat equation is equivalent 
to the Douglas-Gunn scheme and to a scheme given earlier and independently 
by Brian [4] and Douglas [11]. The form that Douglas and Gunn recommend 
for machine computation (see (2.7) of [12]) is not in the delta form, 
but it can easily be rewritten in delta form. The delta form generally 
leads to the most efficient computational algorithm. 

To simplify the comparison of the Douglas-Gunn scheme and the 
formulation of this paper for linear multi step methods, we consider two- 
step methods. The role of the operator Au° defined oy (7.12) is taken 
in the Douglas-Gunn formulation ([12], page 441) by 


(10.2a) 


/ n+1 n+1 

v u - u * . 


where 

(10.2b) u r ; fl - V n + ^u 0 ' 1 
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and 


(10.2c) 


*0 + *1 “ 1 


Here u 


n+1 


is a prediction of u 


n+1 


based on u and u 


n-1 


The choice 


of is ambiguous unless some condition in addition to (10.2c) is 

imposed. The interested reader should refer to the discussion by 
Douglas and Gunn. 

In the fe rnmlation of this paper, the variable to be determined 
from the scheme is E ^p(E)u° rather than u n+ *. It is obvious from (2.5) 
that p(E)u° must be an approximation to At 3u/9t for any consistent 
scheme. Choosing E *p(E)u n as the unknown variable is equivalent to 
imposing the condition 


(10.3) 


n+1 


n+1 


At 


1 + C 3t 


“ + 0(At 2 ) 


in the Douglas-Gunn method. The normalization factor (1 + £) ^ results 
from the parameterization chosen for the two-step method, as can be seen 
by comparing (7.12) and (10.2). By imposing (10.3) and using (10.2b), 
we obtain (10.2c) and, in addition, 


(10.4) 


Hence 


(10.5) 


1 + 2£j 

i + e 


i + f. 
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Consequently, if one epplies a linear two-step scheme [with the param- 
eterization (7.1)] to the heat equation (9.1) and uses the Douglas-Gunn 
method with (♦q*^) determined by (10.5), then the resulting algorithm 
will be equivalent to (8.3). Since condition (10.3) is not part of the 
Douglas-Gunn formulation, the scheme (8.3) will not, in general, 
coincide with the Douglas-Gunn method. It is obvious that if the 
constants (♦q*^) are determined by (10.5), then will depend on 

the particular time differencing method chosen. 

Ue should mention that since Douglas and Gunn did not include 
first-order hyperbolic equations in their formulation, we are demanding 
a more stringent stability condition (namely A-stability) than one would 
require if the eigenvalue [see (4.3b)] were pure real. In the latter 
case, A - A R £ 0 and A^-stability is sufficient for unconditional 
stability. A^-stability means that the region of stability (see 
section 2) contains only the interval (-« ,0] rather than the entire 
left-half plane as required for A-stability. (See also last paragraph 
of section 9.) 

11. Concluding remarks and summary. 

In this paper we have combined A-stable LMMs and approximate 
factorization to construct unconditionally stable ADI schemes for partial 
differential equations with both convection (hyperbolic) and diffusion 
(parabolic) terms. Linear stability analysis for multilevel partial 
difference equations is usually very difficult. The stability of a 
particular scheme is determined by the location of the roots of the 
characteristic polynomial relative to the unit circle in the complex 
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plane. There are testa such as the Schur-Cohn criterion to determine 
when the roots of the characteristic polynomial have modulus less than 
or equal to unity (see, e.g., [21,22]). However, it is often difficult 
in practice to apply thesf. tests because of the complicated nature of 
the coefficients (generally complex) of the characteristic polynomial. 

(See appendix B for the complexity involved for a simple model problem.) 

We have circumvented this difficulty by constructing a class of ADI 
schemes such that a sufficient condition for unconditional stability 
(A-stability) is maintained in the step by step development of the 
scheme. 

The approach is summarized as follows. An A-stable LMM is chosen 
as the basic time differencing scheme. Next, one discretizes in time 
but not in space and applies the LMM method to a (model) linear PDE (4.1). 
The resulting (space continuous) scheme retains the A-stable property 
since the requirement that the real part of the eigenvalue (4.3b) be 
nonpositive is the parabolicity condition for the PDE. The ensuing 
approximations in the construction of an ADI scheme are such that the 
. real part of the "eigenvalue" remains nonpositive. The implicit operator 

I to be inverted is constructed so chat the unknown variable to be 

determined is p(E)u n . This ensures that the approximate factorization 
does not upset either the temporal accuracy (second-order) or the 
stability of the scheme [by lemma (5.7)]. Finally, one notes that 
(central) spatial discretizations do not alter the essential property 
for A-stability, i.e., the real part of the eigenvalue is nonpositive. 

Since our main emphasis in this paper is linear stability theory, 

] we have considered only model equations. In a companion paper [3], 
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we apply the method outlined In section 3 to derive an ADI algorithm 
for a mixed hyperbolic-parabolic system of nonlinear equations where 
the time differencing is the class of A-stable linear two-step methods. 
Earlier references on the development of noniterative ADI schemes for 
nonlinear systems of partial differential equations are listed in the 
introduction. 


% 
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Appendix A. Stability analysis when model equation includes a mixed 
spatial derivative. 

In the stability analysis of section 5 we assumed that the mixed 
derivative term bu^ the model equation (4.1) was zero. In this 
appendix we show that if the mixed derivative is treated explicitly by 
a first-order method, then the resulting ADI scheme remains uncondition- 
ally stable. An explicit treatment of the mixed derivative means that 

bu does not appear on the left-hand side of the factored scheme (3.3) 
xy 

but does appear on the right-hand side. Let the mixed derivative be 
appended to (3.3) as follows: 


(1 - niAt L )(1 - wAt L )p(E)u = 
x y 


At[o(E) - wp (E) ] IL + L + b 

V x y 3x3y 


3i\n 
x3y/ U ’ 


where L and L are defined by (5.1). 
x y 

By following the same stability analysis that led to characteristic 
equation (5.5a), we obtain 




(A2a) 


p(C) - A At o(;) = 0 , 


where A,„ is now defined bv 
Ar 


(A2b) 


A j + _ hk i K 2 


1 + uAu A^A ” W ‘' it 


and A^ and A^ are defined by (5.5b) and b is the coefficient of 

u in the partial differential equation. A sufficient condition for 

xv 
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unconditional stability is that the Re £ 0. Since has the 

AF AF 


°i + ie i 


AF a 2 + iB 2 * 


Re X A p £ 0 if and only if 


a l a 2 + S 0 * 


By a direct calculation we find 


+ B^B 2 “ - ^ K i 2 + ' 3 ic 1 k 2 + cic 2^^ ac ^^ “ bi|» + 1^ 

2 2 2 / 2. . 2\ 

‘ “ ‘l *2 \ cc 1 - bc l c 2 * ac 2 ) • 


where <J» ■ w At k^< 2 . quadratic forms 


2 2 
a<^ + bK^< 2 + ck 2 


2 . , 2 
cc t - bc^ + ac ? 


are nonnegative if and only if a,c > 0 and b" < 4ac, which iV the 
parabolicity condition for the l’DE. (If equality holds the equation is 
hyperbolic.) Likewise, the quadratic 



F»* 


... till 1 1 1 l I_l 1 1 1 i l I I Pa Li. 


2 

is nonnegative if and only if b £ 4ac. Consequently, Re A^p £ 0 and 
the ADI scheme with the mixed derivative term treated explicitly is 
unconditionally stable. The above analysis is for the spatially 
continuous case. However, in the spatially discrete case, A^ will also 
have the same form as (A2) and thus the spatially discrete ADI scheme 
will be unconditionally stable. 

It is very interesting to observe that if one treats the mixed 

derivative term explicitly but does not use a factored scheme, the 
2 2 

product term u> At A^A^ would be missing from the denominator of (A2b). 

In this case one can always pick wave numbers 80 t * lat t * ie 

denominator is negative, and then Re A > 0. Consequently, the 

Ar 

unfactored scheme is not unconditionally stable. 

Finally, we should point out that if the cross derivative term is 
treated explicitly as above, then the resulting scheme is first-order 
accurate in time for the mixed derivative. (See last paragraph of 
section 5.) 
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Appendix B. Stability analysis for an alternative formulation. 


In this appendix we carry out a linear stability analysis for the 
factored scheme (9.4). We assume a solution of the form (5.2) and find 
that the Fourier coefficient satisfies 


(Bl) (1 + ic^d + ic^Av" - - | (ic x + ic 2 )v n 


- | lie, + i^lAv"" 1 + 


where we have defined 


l' At - 0 At 

C 1 ~ 1 + f. Vl ’ c 2 “ 1 + t; V 2 


and used the Identity 


„ n * n-i 

Vv * Av 


From (Bl) it follows that the amplification factor defined by (5.4) 
satisfies the quadratic equation 


a 2 C + V* + a o " 0 • 


whe ri- 


Cl + i-' j ) ( 1 + ic , ) . 


-(! + iCjXl + i" 2 ) +(~ + + ~ 


(i ^l + i ^2 ) + 7 + 
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The characteristic equation (B4) with coefficients (B5) cannot be p«it 
in the form (5.3a) and, consequently, we must use some other method of 
determining when the modulus of the roots of (B4) are bounded by unity. 
Here we use the Schur-Cohn criterion 121] as formulated by Miller [ 22 j . 

The polynomial (B4) is a von Neumann polynomial [22], that is, 

| c | SI, if and only if either 


(B6a) 


a„a„ 
0 0 


a 2 a 2 


0 


and 


(B6b) . 


A 2 * <Vo 


- a 2 J 2 > - (.i 2 a ( 


Vl )(;, 2 a l 


Vi> 


or 

(B6c) A 1 " A 2 " 0 ’ 4a 2 a 2 " a i a l * 0 * 

Substitution of the coefficients (B5) into (B6a,b) yields, after some 
algebraic manipulations. 
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*2 " “e — ~] t* 1 + Ot<26 + 2* - 1) - £ (26 + 24> 1)J 


2 

+ (20 + 2 $ - 1 )^ 2 c 2 2 + (“-^) (26 


2 $ - l)^ 2 


+ 2 l - (26 + 2<() - 1) + 


(4-1) (28 - 2* - l)]v 


(rr*) (2e - 2* - d : 2 2 


where £ ■ £/(l + £)• For the class of all two-step methods (9.2) that 

are at least second-order accurate, the parameters (0,£,4>) are related 

by (7.3). For this class of methods we want to determine the parameter 

space (6,0 for which the factored scheme (9.4) is unconditionally stable. 

First, we consider the conditions under which as given by (B8) 

is nonnegative. Using (7.3) we find that the expression enclosed in 

the first pair of brackets within the braces of (B8) is zero. The 
- 2-2 

coefficient of c^ c^ is nonnegative if and only if 


C > 0 


The remaining tents within the braces constitute a quadratic form 


(BIO) 


- 2 - 2 

Ac, + Bc.c + Cc , 

1 i 2 2 


which is nonnegative if and only if 


(Blla,b,c) 


A > 0 , C > 0 , B < 4AC 


The first two inequalities (Blla.b) require 


1 





i_j ^ 1 1 *1 l ! l — L a— 1 — L - i 


1 . 


(B12a) 


C i 29 - 1 . 


Necessary and sufficient conditions for the third inequality (Bile) 
are 

(B13a) 0 ii (1 + t)[(l + 20 - V (1 + 20(1 + o]/€ , 

(B13b) 0 Ml + O [(1 + 20 + Vd + 2 0(1 + Oj/5 

for £ > 0. 

Formula (37) for A^ can be rewritten as the sum of two quadratic forms: 


-_ d i 


d 2 c l c 2 


(V 2 ) 




+ Vl'2 


d 3 J 2 2 


]• 


where 



Necessary and sufficient conditions to satisfy < 0 are less 
restrictive than conditions (B9) and (Bl3a,b), as the interested reader 
can verify. Hence, the region of unconditional stability is defined by 
the inequalities (B9) and (B13) and is indicated by the shaded region 
of figure 3 of the text. 
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